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® A method for processing seismic data. 

@ A method for processing seismic data wliich include performing normal moveout velocity analysis in such a 
manner as to include the effects of the offset dependence of reflection amplitude and simultaneously determin- 
ing from the seismic data the moveout velocity and the variation of reflection amplitude with offset. 
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A METHOD FOR PROCESSING SEISMIC DATA 



The present invention relates to exploration seismic reflection surveying and more particularly relates to 
the proce.ssing of exploration seismic reflection data to enliance information in seismic signals reflected 
from contrasts or differences in elastic constants or densities In the subsurface of the earth. 

The method of the present invention which is described herein is generally discussed In terms of 

5 compressional-wave (PP) seismic data acquisition and processing, which is the most common form of 
seismic data used in exploration seismology. However, it should be understood that the method is equally 
applicable to shear-wave seismic data. 

Conventional land or marine seismic acquisition techniques involve the use of an appropriate source to 
generate seismic energy and a set of receivers, spread out along or near the surface of the earth on land or 

10 at or near the water surface or water bottom in a water covered area, to detect any reflected seismic signals 
due to seismic energy striking subsurface geologic boundaries. These signals are recorded as a function of 
time and subsequent processing of these time varying signals, i.e. seismic "traces" or seismic data, Is 
designed to reconstruct an appropriate Image of the geologic boundaries of the subsurface and to obtain 
information about the subsurface materials. In simplistic terms, this conventional process has a seismic 

75 wave, from a source of seismic energy, travelling down into the earth, reflecting from a particular geologic 
interface (i.e. a change or contrast in elastic constants and/or densities), and returning to the surface, where 
It may be detected by an appropriate receiver. 

If the seismic-wave velocity Is known as a function of depth and lateral position, and if the position and 
dip of a plane geologic Interface are known, the time for the wave to travel down to that particular reflecting 

20 interface and reflect back to the surface can be computed for any source and receiver locations. This two- 
way travel time is usually described by a function tp<.Z). where Z is the depth to the reflecting interface 
(contrast in elastic constants or density) and X is the horizontal distance (offset) between source and 
receiver. 

If the elastic constants and densities of the materials above and below a planar reflecting interface are 

25 known, then the reflection coefficient may be computed. This reflection coefficient is the ratio of reflected 
amplitude to incident amplitude and will depend on the angle of incidence at the reflecting interface. The 
angle of incidence, 0. is the angle between the ray normal to the incident downgoing wave-front and a line 
normal to the interface; as is well known, the incident and reflected rays will be in a plane normal to the 
interface. This angle of incidence increases with Increasing offset X. The reflection coefficient for a 

30 compressional wave from a particular interface will be designated by the function Rp (e), where the angle 8 
may be related to the offset distance X and depth of reflector Z by raytracing if the compressionalwave 
velocity at all points In the earth Is known; this velocity information, or a reasonable approximation thereto, 
is referred to as a "velocity model". For a given reflector, the reflection angle, 8, and offset, X, are 
geometrically related, so any discussions herein in terms of dependence upon offset (offset dependence) is 

35 equivalent to dependence upon reflection angle (angular dependence). The angular (or offset) dependence 
of reflection amplitude may be computed exactly for a point source and plane reflector, however In most 
practical cases it may be approximated adequately by plane-wave reflection coefficients (reflection coeffi- 
cients for an incident plane wave) which are easily calculated using expressions derived from the results of 
Zoeppritz (see for example. Bull Seismoi. Soc. Am., Vol. 66, 1976, pp. 1881-1885, Young. G. B. and Braile. 

40 L.W.). For a compressional-wave reflection from a planar interface between two media having a small 
contrast (i.e., with the medium containing the incident and reflected waves having a compressional velocity 
Vp . a shear velocity Vs , and a density p, and the other medium having a compressional velocity of Vp + 
dVp , a shear velocity Vs + dVs , and a density p + dp, and where dVp A/p , dVsA/si and dp/p are small 
compared to one), the offset (or reflection angle) dependence of reflection amplitude may be described for 

45 angles of incidence less than the critical angle by an expansion of the form, 
Rp(0) « Rp(0) + Ksin2(o) + Lsin*(e)+ (1) 

For the discussion herein, the angles of Incidence are limited to angles such that the terms of the order of 
sin* (0) and higher are negligible. In Equation (1), Rp (0) Is the normal incidence {8-0) reflection 
coefficient; Rp (0) depends only on the densities and compressional velocities of the two media. K Is a 
50 constant, which also depends on the elastic properties and densities of the media. The relationship of K to 
the elastic properties and densities may be expressed in a number of ways. One particulariy simple 
expression which relates K to the contrasts in shear velocities, compressional velocities, and densities is K 
= Ra-4(VsA^p)2(2R/j + Rp, (2) 
where 

Ra=dVp/(2Vp + dVp). -(2a) 

2 
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Rj3 = d Vs/(2Vs + dVs). and (2b) 
Rp = dp/(2p + dp) 
Also, in terms c 
given exactly by, 



Ro = dp/(2p + dp). (2c) ^ . ^ 

Also, in terms of these same coefficients, the normal incidence or zero offset reflection coefficients are 



TO 



^P^°^= "T^' and (2d) 



R (0)= P P . (2e) 



S l+RpRp 



For sufficiently small values of Ra. R/3. and Bp . equations 2d and 2e may be approximated as, Rp(0) = R« 
^ + Rpand (2f) 

Thus m^urement of ttie normal incidence compressional-wave reflection coefficient, Rp (0). gives 
information about the densities and compresslonal velocities, v^hile measurement of the offset dependence 
constant K can provide information about the densities and shear velocities of the media. The formulas 
given above are for small contrasts in the elastic properties and densities above and below the planar 

25 interface- of course, more general theoretical relations may be used. Similar relationships are well known for 
the offset dependence of shear-wave reflection coefficients, although the particular form for such shear- 
wave equations that are analogous to equations 2-2g Is quite different. 

There are a number of geologic questions important to exploration for hydrocarbons which can be 
answered by acquiring a knowledge of both the compresslonal- and shear-wave properties of the subsur- 

^ face materials. For instance, these materials are generally porous with various fluids filling the pore space. 
The compresslonal velocity of a seismic wave in such media depends strongly on the rock matnx 
properties as well as those of the pore fluid. On the other hand, shear-wave velocities depend strongly on 
the rock matrix but only slightly on the pore fluid. Thus, detailed study of both the compresslonal and shear 
properties of the media provides an opportunity to characterize any changes in seismic response as being 

35 due to changes in fluid content (e.g. from brine to oil. or oil to gas) or changes in the rock matrix (e.g. from 
sandstone to shale or a change in porosity). The ratio of Vp to Vg is often a useful diagnostic feature of such 
changes It should be noted that, even without lateral variation. In many cases the recognition of fluid 
content or rock type may be possible with an accurate knowledge of the compresslonal and shear 
properties at a single location. Distinguishing between fluid effects and lithology effects, and detecting 
different porosity and lithology types are of vital exploration interest and the desire to make such 
distinctions has engendered significant effort in the measurement and Interpretation of shear properties in 
addition to the information concerning compresslonal properties traditionally inferred from conventional PP 
reflection prospecting. 

It is generally the objective of seismic exploration to generate seismic energy, make measurements of 
« the reflection amplitude of this energy at various offsets and for various times, and then, by employing 
various processing steps on this seismic data, to deduce the geometry as well as some of the elastic 
properties and densities of the materials of the earth through which the seismic energy has propagated and 

from which it has been reflected. ^ ...... , •.. 

One such processing method which resulte in some knowledge of a velocity model (that is the velocity 

^ of seismic waves as a function of Z) Is conventional moveout velocity analysis. This velocity model may 
then be used for dynamic correction followed by stacking the corrected seismic data, as described below. 

A seismic reflection from an interface will arrive at a receiver after a two-way travel time, denoted and 
used herein as. t(X). where X is defined and used herein as the distance between the source and the 
receiver or "offset" distance. This "moveout time" t(X) may be used to 'dynamically correct" seismic data 

^ acquired at an offset distance X so that a reflection Is adjusted in time to appear as if it had been acquired 
at zero offset, i.e. X=0. Conventional "stacking" is accomplished by summing or averaging such dynam- 
ically conrected data. 

For seismic data where the downward and upward propagating wave types are the same wave type. 
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either compressional or shear, the square of this moveout time has a well known expansion about X=0: 

= to2 + 32X2 + 32X2 +a4X* + (3) 

where the coefficients aa, a*, ae, etc. ... are known functions of the propagation velocity at each depth and 
to is the X = 0 (zero offset) two-way travel time. The propagation velocity at a particular depth Is known as 

5 the "interval" velocity at that depth. Experience with compressional-wave data has indicated that in most 
areas the P-wave vertical velocity distributions are such that, within the offset range (X < Z) and frequency 
range {5-100Hz) normally considered in exploration seismology, terms of higher power than X2 are small 
and for many purposes may be Ignored. Ignoring the higher order terms is equivalent to approximating the 
stratified earth with a single horizontal layer having an "effective" velocity, Vg . For this simplified case 

TO Equation 3 reduces to the "normal hyperbolic moveout" expression t2p<) = to2 + X^/Ve^. (4) 
where 

2 . ' J '^'^ (4a) 
1 • 



20 



25 



30 



and 

to = 2 / [dZA/(Z)]. {4b) 

For the simplistic discussions hereinabove, the earth model is assumed to be one for which velocities 
may vary with depth, but not laterally; however, the method of the present Invention may be applied in more 
general cases. For example, this model and analysis is also appropriate for reflections from dipping 
reflectors as long as the data is gathered by common midpoint (CMP) between source and receiver pairs 
(as more fully described below), and the degree of dip is not severe, and the reflectors are laterally 
continuous, if these assumptions are not met, other conventional processing steps, such as migration, may 
be required, although they are not described herein. For a dipping reflector, the interpretation of the 
coefficient aa of equation 3 as the square of the reciprocal of effective velocity must be modified to include 
the cosine of the dip angle. For this case equation 4 becomes: 

X^X)=t^ + 5 , (4-) 

3S 

where o is the dip angle. 

Conventional processing of compressional-wave data uses data collected with many sources and many 
receivers and gathers the data by the common midpoint (CMP) technique, as illustrated in Rg. lA of the 
accompanying drawings, which will be discussed in more detail hereinafter. Traces with a common 

40 "midpoint" between the source and receiver are collected or gathered at a surface gather point. For 
example, in Fig. lA, Si and Ri are the source and receiver pair for the first trace and have a midpoint at 
the surface point (0). Rg. 1B of the accompanying drawings, which will be discussed in more detail 
hereinafter depicts the conresponding hyperbolic moveout of such data (where the numbers used cor- 
respond to the subscripts used in Fig. 1A) and Fig. 1C of the accompanying drawings, which will be 

45 discussed in more detail hereinafter depicts the corresponding variation of reflection coefficient with offset 
for such a case. That is. after the data have been acquired for a number of sequential source positions the 
traces for various source-receiver combinations are sorted or gathered into different midpoint groups which 
have the same or "common" surface location of the "midpoint" between the source and receiver positions 
(as depicted in Rg. 1A). This sorted or gathered midpoint data Is then analyzed or processed (via equations 

50 4) to determine effective velocities Ve , also called normal moveout velocities, for reflections from various 
depths, i.e. various values of to. One method is to determine for each to a value of Ve which provides 
dynamic con-ections which maximize the resultant amplitudes of the "stacked" data in a time gate around 
to. Many methods for moveout velocity analysis and dynamic correction of compressional (PP) and shear 
(SS) seismic data have been in use since the late 1960*s, as described for instance by Tanner and Koehler 

55 ("Velocity Spectra -Digital Computer Derivation And Applications Of Velocity Functions", Geophysics 34, 
pp. 859-881 (1969) Tanner, M.T. and Koehler. F.). 

The original basis for CMP processing is the fact that each trace in a gather images (or consists of 
reflections from) approximately the same subsurface points, and, when properly adjusted for differing path 
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lengths, the set of corrected traces may be averaged to give an enhanced picture of the reflection response 
of the earth below that CMP surface location by emphasizing taie primary reflections and discriminating 
against multiple reflections and other undesirable noise. It is usually assumed that the "stacked" trace 
represents the normal incidence (zero-offset) response of the earth. While this procedure has been very 
6 effective in improving signal-to-noise ratios for seismic data in many areas, it ignores the fact that reflection 
amplitudes vary as a function of offset and that the stacked trace is not equivalent to a normal incidence 
trace. 

More recently, methods have been described for measuring and interpreting the variation with offset of 
the reflection amplitude from a given subsurface interface. Techniques which are taught for example in U.S. 

10 Patent Specifications Nos. 4.562.558; 4.573.148; 4,570.246; 4.316.267; 4.316.268 and 4.534.019 explicitly 
describe methods for measuring and interpreting amplitude variation with offset Since reflections appear at 
increasing times at increasing offsets, all these newer methods also require a technique for identifying the 
time at which a given reflection appears at each offset. Various conventional methods have been used for 
identifying these times in the above-noted patents, as well as in other publications on the subject. One 

15 method simply proposes following a particular trace feature, say a maximum amplitude as it moves out with 
offset. A more common method utilizes formulas like equation (4) with an effective velocity determined by 
conventional means (i.e. ignoring amplitude versus offset) as described in the references noted 
hereinabove. 

Now it has been discovered that in many cases neither the conventional technique of event moveout 

20 determination by conventional velocity analysis, or the conventional technique of following the peak 
amplitude (or other feature) of an event as a function of offset is of sufficient accuracy for use in studying 
reflection coefficient variation with offset, specifically because these techniques do not properiy account for 
the effects on moveout determination of amplitude variation with offset. Examples, where these conventional 
techniques fail, include situations where the reflection amplitude varies with offset in such an extreme 

25 manner that \he reflection amplitude changes sign, going from positive values to zero and then to negative 
values over the ranges of offsets being considered. Another case is described later herein, where the 
varying offset dependence of closely spaced reflectors can have the appearance of a single reflector 
convolved with a wavelet which changes with offset in a spurious- manner. Thus, a method for determining 
the effective velocity of seismic events more accurately by specifically including the effects of amplitude 

30 variations with offset is needed. Moreover, such a method could simultaneously determine a more accurate 
measure of amplitude variation with offset and normal incidence reflection amplitude. 

These and other limitations and disadvantages of the prior art are overcome by the present invention, 
however, and an improved method is provided for determining from seismic data normal moveout velocities 
and associated dynamic corrections that preserve amplitude versus offset information to provide better 

35 estimates of effective velocities, the nonmal incidence reflection amplitudes, and the offset dependence of 
the reflection amplitudes. 

It is an object of the present invention to provide a method for performing normal moveout corrections 
which allows for determination from seismic data of true reflection amplitudes at zero-offset or normal 
incidence. 

40 It is another object of the present invention to provide a method for performing normal moveout 
corrections which allows for determination from seismic data of true reflection amplitude dependence on 
offset. 

It is another object of the present invention to provide a method for performing normal moveout 
con-ections which allows for determination from seismic data of true reflection amplitudes at zero-offset and 
45 true amplitude dependence on offset. 

It Is another object of the present invention to provide methods for performing moveout velocity analysis 
that include consideration of appropriate amplitude variation with offset 

It is another object of the present invention to provide methods for more accurately determining 
effective velocities, from which more accurate interval velocities may be derived. 
60 It is another object of the present invention to provide methods for constructing more accurate and 
more easily interpreted seismic sections. 

It is another object of the present invention to accurately determine the normal moveout velocity, normal 
incidence reflection coefficients, and the offset dependence of the reflection coefficients, which may then be 
used to provide information about the elastic constants and densities of the subsurface materials which 
55 have caused the reflection of seismic signals. 

It is another object of the present invention to provide a method for processing seismic data, 
comprising, simultaneously determining the moveout velocity, normal incidence reflection amplitude, and 
offset dependence of reflection amplitude from said data. 
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It is another object of the present invention to provide a method for processing seismic data, 
comprising, calculating from a known reference velocity function effective reflection angles for various two- 
way travel times and offsets of said data, detemiining from said data and effective reflection angles for a 
trial moveout velocity function, fitting coefficients for a preselected offset dependence function which 

5 expressly accounts for offset dependence of reflection amplitudes, repeating said determining of said fitting 
coefficients for a plurality of said trial moveout velocity functions, and selecting for various normal incidence 
two-way travel times appropriate moveout velocities and fitting coefficients from said data. 

It is another object of the present invention to provide a method for determining properties of sub- 
surface materials, comprising, generating seismic energy on or near the surface of the earth, measuring on 

10 or near the surface seismic signals resulting from said seismic energy reflected from interfaces, determining 
from said measured seismic signals optimized moveout velocities, normal incidence reflection amplitudes, 
and angular dependencies of reflection amplitudes, and determining elastic properties and densities of said 
subsurface materials from said moveout velocities, normal incidence reflection amplitudes, and offset 
dependencies of reflection amplitudes, 

75 Therefore, the invention provides a method for processing seismic data to enhance information in 
seismic signals reflected from contrasts or differences in elastic constants or densities in the subsurface of 
the earth, comprising the steps of: 

- generating seismic energy on or near the surface of the earth; 

• measuring on or near said surface seismic signals resulting from said seismic energy reflected from 
20 interfaces, and recording said signals to obtain seismic data characterized by the step of specifically 
including the effects of amplitude variations with offset. 

In an advantageous embodiment of the present Invention seismic data that has been acquired 
employing a plurality of l<nown source and receiver locations is processed to provide simultaneously 
determined and self-consistent representations of a) zero-offset reflection amplitudes, b) reflection amplitude 
25 dependence on offset, and c) effective velocities, all determined from the seismic data. The method of the 
present invention employs amplitude versus offset Infoanation during such processing, to provide more 
correct velocity and dynamic correction information than Is available from conventional techniques and to 
provide simultaneous quantitative estimates of amplitude variation with offset. 

These and other objects and advantages of the present invention will bec9me apparent from the 
30 following detailed description, wherein reference Is made to tiie Rgures In tiie accornpanying drawings. 

Rgs. 1A, 1B, and 1C depict conventional seismic common midpoint acquisition geometry, the hyper- 
bolic moveout of reflection data obtained, and the variation of reflection coefficient with offset, respectively; 

- Figs. 2A-2D depict a model of a layered earth, associated reflection coefficients at various offsets, and two 
synthetic seismograms derived from this model; 

35 - Rg. 3 depicts a simplified flow chart of the processing steps of the method of the present invention; 

- Rg. 4 depicts synthetic common midpoint reflection data dynamically corrected with con'ect velocities; 

- Rgs. 5A and 5B depict, for the data of Rg. 4, simplified representations of quality of fit coefficients of 
conventional velocity analysis and velocity analysis of the present invention, respectively, plotted as a 
function of t and estimated Ve; 

40 - Rgs. 6A and 6B depict, for the data of Rg. 4, quality of fit coefficients recalculated and plotted as a 
function of t and the difference between estimated Vq and true effective velocities from the model; 

- Rgs. 7A and 78 depict panels of A(t) and B{t) traces, respectively, corresponding to the individual velocity 
functions used to generate the quality of tit coefficients In Rg. 6B; 

- Rg. 8 depicts enlarged 6(t) traces and corresponding model determined coefficients of sln^e, and 

45 - Rgs. 9A and 9B depict, respectively, A(t) and A(t)-B(t) panels for real seismic data recorded along a 
particular line of common midpoints. 

As already indicated earlier the present invention provides a method for processing conventionally 
acquired seismic field data. The method includes normal moveout velocity analysis, including appropriate 
offset dependence, and dynamic conrection of tiie data. The velocity analysis uses offset dependence 

50 parameters determined simultaneously from the data. Additional processing steps may be performed as 
necessary or desired before or after any of these steps. For example, the data may also be corrected for 
statics or othenArise manipulated to enhance signals. 

The metiiod of the present invention which is described herein are generally discussed in terms of 
compressional-wave (PP) seismic data acquisition and processing, which is tfie most common form of 

55 seismic data used in exploration seismology. It should be cleariy understood that these methods are equally 
applicable to shear-wave seismic data acquired using an appropriate shear source and appropriate shear 
motion detectors. Other than the particular acquisition techniques involved and the techniques for final 
interpretation of seismic attributes, tiie only changes necessary for the application to shear data would be to 
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replace the discussions of the compressional-wave moveout and interval velocities with appropriate shear- 
wave moveout and interval velocities. While the measurement of normal moveout velocities, nomnal 
incidence reflection amplitudes, and the offset dependence of the reflection amplitudes, may be accom- 
plished in the same manner for each of these modes of seismic propagation, the relation of these 
5 measurements to the elastic properties and densities of the medium will not be the same. The reflection 
coefficients for these various modes have well known but quite different dependences on the elastic 
constants and densities of the reflecting materials. 

The following phrases or expressions used herein, are defined as follows: "Amplitude versus offset 
(AVO)" means the phenomenon of seismic signal amplitude varying with source-receiver offset; such 
TO variation is caused by angular dependence of reflection coefficients, wavefront divergence, acquisition 
techniques, absorption, etc.; "Offset dependence of reflection amplitude" means that part of AVO due 
specifically to reflection coefficient variation with offset; and "Angular dependence of reflection amplitude" 
or "Angular dependence of reflection coefficients" means the variation of amplitude or reflection coefficient 
with angle of incidence, an effect which may be deduced from AVO by various processing and measure- 
15 ment steps along with raytracing to relate reflection angle to source-receiver offset. The angular depen- 
dence of the reflection coefficient is ultimately required for determining elastic constants and densities of 
the materials Involved in reflections. All of these expressions are related and in many statements referring to 
a common property of the three terms they may be used interchangeably. Further, the following terms are 
also used as follows herein: "Normal incidence" means zero offset (X=0) for PP. and SS data; "PP" means 
20 a downgoing and upcoming compressional wave, with particle motion parallel to direction of propagation; 
and "SS" means a downgoing and upcoming shear wave, with particle motion perpendicular to the direction 
of propagation. Also, note that for notational convenience the two-way travel time at zero offset, defined 
hereinbefore as to , is hereinafter denoted as t. 

Referring now to Figs. 2A-D. it may be seen that conventional moveout velocity analysis may result in 
25 erroneous velocity estimates. Improper dynamic correction, and incorrect amplitude interpretation of 
particular seismic events. Fig. 2A represents a simple model of two 16 metres Intervals of a sand layer, one 
gas filled (indicated by A) and one brine filled, embedded in shale (indicated by C and D respectively) at a 
depth of 3300 metres. Realistic velocities and densities were used in Zoeppritz's equations to determine 
reflection coefficients for plane waves at the angles of incidence at the reflectors corresponding to various 
30 source-receiver separations, X. The reflections from the three interfaces were found to be separated by two- 
way travel times on the order of 10 msec for X = 0; these reflection tweway travel time separations 
gradually decrease as X increases. Rg. 2B depicts the offset dependence of the reflection coefficient for 
each individual reflector or interface. The horizontal axis represents X offset, whereas the vertical axis 
represents amplitude. The plotting symbols correspond to the symbols identifying the interfaces in Rg. 2A. 
35 The displays in Rgs. 20 and 2D show a number of synthetic seismogram traces, computed for various 
source-receiver offset distances. The traces are plotted vertically with time (in milliseconds) Increasing 
downward along the vertical axis; they are spaced along the horizontal axis at distances proportional to the 
offset. The synthetic selsmograms represent the expected seismometer responses to earth motion with 
motion in the direction defined as positive representing a trace deflection to the right, A given motion of the 
40 earth would produce the same amount of deflection in each synthetic trace. 

Each trace in the synthetic selsmograms depicted in Rgs. 20 and 2D was constructed by convolving a 
waveform with a reflectivity function made up of three impulses (one for each interface). The amplitudes of 
the impulses are the plane-wave reflection coefficients for the offset indicated for the trace, and the 
uncorrected times at which they were posted are the travel times of a compressional wave from the source 
46 to the reflector and then to the receiver at the offset indicated. Note that for the waveform chosen, which 
represents a realistic seismic frequency range, the impulses are not resolved as separate events. 
Rg. 20 depicts the synthetic traces dynamically corrected with ttie conrect moveout. or effective velocity, 
associated witii tiie reflection from each interface at tills depth. The connect velocity value is computed from 
tiie model. Rg. 2D depicts the synthetic traces dynamically corrected with a moveout velocity function 
50 chosen by conventional moveout velocity analysis to maximize the correlation of the dynamically corrected 
traces. Although tiie correct effective velocity for the middle reflector is 7794 (fps) (266 metre per second), 
tiie conventional moveout velocity analysis provides a value of 7700 (fps) (233 metre per second). Thus, tiie 
conventional moveout velocity analysis provides an erroneous estimate of the taie effective velocity. 

Note that a measurement of AVO made at a time of 2600 milliseconds on Fig. 2D would give a very 
55 different result than a similar measurement made on Rg. 20 where correct effective velocities were used. 
That is, at an offset. X. of zero on botii figures, the same maximum right deflection occurs at the timing line 
while at an offset, X, of 10.000 ttie correct amplitude at the timing line in Rg. 20 is much smaller than tiie 
amplitude seen In Rg. 20. 
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10 



Before giving a detailed description of the present invention's method of velocity analysis, a detailed 
description of conventional velocity analysis is provided. Conventional velocity analysis involves a dynamic 
correction process using an assumed trial moveout velocity function Ve (t), followed by a measurement of 
the deviation of the corrected data from an assumed Ideal situation in which properly corrected reflections 
have the same amplitude and shape for all offsets (which is the same as maximizing the amplitude of the 
stacked trace or maximizing the correlation of ail the dynamically corrected traces). The measurement of 
the deviation may be done on a (time) sample by (time) sample basis, although normally some averaging 
over adjacent times is applied. For a particular assumed Ve (t), a "quality of fit coefficient". E[Ve (t)], which 
depends on the squared differences between the trial corrected data and the ideal, may be computed and 
then plotted versus velocity and two-way travel time. Note that E[Ve (t)] is given by, 



IS 



E(V^(t)l = 1 - 



I min I (C'(t»)-S(f ,X)) 
t' C * 



(5) 



the summation over X is over all offsets, the summation of t Is over samples in the time window between t- 
At and t + At . S(t.X) is a seismic panel dynamically corrected for the trial velocity (t), and "min" 
20 indicates that C'(t') is chosen to minimize the sum of the squared error differences for each time t . (Of 
course, it is also possible to only employ the bracketed term in equation 5 and have quality of fit 
coefficients that are minimized, i.e. most nearly approach zero.) The step of computing and plotting the 
quality of fit coefficient E[Ve (t)] is repeated for various Ve (t) functions to determine what velocities provide 
th e l ar g est qual i t y-o f - &t - GQ6ff i cient ^at■a .g i ver^ tlme r t -- 



2^ For the hereinabove described case, the sum of squared errors Is minimized for a particular t when C 
= (I 'Nx ) § S(t.X). where Nx is the number of offsets with valid data, and equation 5, reduces to 



1 I llsit\x)V 
:[vjt)i = t' \ ^ 



EI . 

30 « ^ (5a) 



Equation 5a is also an expression for the relative power of the stacked trace and also the correlation of all 
the dynamically corrected traces in a time window about t. Thus, maximizing the quality of fit coefficient 
also yields the maximum relative power and maximum correlation. 

In order to include the effect of amplitude variation with offset in the methods of the present invention 
for velocity analysis such amplitude variation must be characterized. As noted hereinbefore, the general 
characteristics of offset dependence of reflection amplitudes have been found to be adequately approxi- 
mated for a given interface as: 
Bp(e) ' Rp(0) + K sin2(9) (6) 

where Rp (6) is the normal incidence (0 = 0) reflection coefficient and K is a constant, and both Rp (0) and 
K depend on the elastic properties and densities of the two media, as described in equations 2. 

Although the particular term used for including the offset dependence of reflection amplitude is 
preferably sin^ (d), other terms may be employed In the methods of the present invention. Alternative 
approximation expressions such as. 
Rp(9) = Rp{0) + k' tan2(0) 
or 

Rp(0) = Rp(0) + k" X2 

may work nearly as well in fitting the offset dependence obsen/ed in most data over a limited offset range. 

The approximation expression employed to describe the offset dependence will have some effect on how 
large an offset range may be used in the processing and how large the reflection contrast may be before 
the approximation becomes inadequate. The choice of offset variable (i.e., source-receiver separation, 
reflection angle, injection/emergence angle, etc.) has littie effect on tiie quality of fit to tiie data, but does 
affect how tine measurement of the offset dependence can be interpreted in terms of tiie elastic properties 
and densities. A con'ect interpretation requires the use of tiie correct reflection angle which can be 
determined from tiie known offset distance X and travel time t If tiie exact interval velocity model is known. 
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Since this velocity model is seldom known exactly, approximate velocities, and therefore approximate 
reflection angles are often determined from the moveout velocity analysis. However, the moveout velocity 
determination and the norma! incidence reflection coefficient determination are both insensitive to the 
approximation of the reflection angle, as well as to the particular term used for including the offset 
5 dependence of the reflection coefficients. , 

The methods of the present invention substitute a function of the form A(t) + B(t) sin2fl{X) for C (t) in 
equation 5, i.e. 

I oia I [A(t') + B(f ) sin^eCX) - S(f ,X)1^ 
E[V ft)) = 1 - A.B 



(7) 

The "min" values of A and B are time dependent fitting coefficients, with A(t) representing a least-squared- 
errors estimate of normal incidence (zero-offset) amplitudes and B(t) representing a least-squared-errors 
measure of the offset dependence of the seismic data for a given Ve (t) function at a time, t. This procedure 

20 Is repeated for various Ve (t) functions. The Ve (t) for which the quality of fit is maximum at a time t is 
chosen as the optimum value and the corresponding A(t) and B(t) are similarly called "optimal" values, as 
described in more detail hereinbelow. 

The summations over offsets in equations 5 and 7 are typically limited to a range corresponding to 
some preselected angle limits or a preselected range of X^. The limits serve several purposes and are 

25 selected to be compatible with these purposes: 1) frequently these limits are chosen to avoid coherent 
noises from arrivals other than primary reflections; 2) they avoid "apparent wavelet stretch" introduced by 
dynamic corrections and due to the convergence of moveout hyperbolas; and 3) they avoid the angles or 
offsets so large that the seismic energy is critically reflected or the approximate form in equation 6 is not 
sufficiently accurate. Although items 1) and 2) above are Important for conventional analysis, all three 

30 purposes are important for the present invention. Range limits are often imposed by "muting" the 
dynamically corrected data, that is, by setting the seismic amplitudes outside the acceptable range to zero. 
The muted data are not counted or included in any of the summations over offset or time in alt of the 
previous and following equations described herein. 

Minimization of the sum in equation 7 with respect to A(t) and B(t) Is a least-squared-ent)r optimization 

35 problem. Stating the problem in terms of a "least-squared-error" optimization allows for convenient 
generalization, for instance, by including various weights based on offset or signal-to-noise ratio. In the 
least-squared-enror method explicitly described herein the weights are taken to be equal. Including more 
general weights is considered within the scope of the present invention. The solution and the resulting 
optimal values may be expressed as simple moments^ of the seismic data, as a function of angle (or 

^ equivalently offset). For equation 5 the optimal value of C may also be written as, 
C' (t) = <S(t.X)»x. (8) 

where the following notation ("< >") is employed for moments of the seismic data 



1 



<S(t,X)sin''(e)>X = - p(t,X) sia'^O). (9) 



SO 

Similarly for equation 7, the optimal values of A and B may be written as: 



55 
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(lOa) 

<S(t,X)>„<sin*e>y - <S(t,X)sin28> <sia2e>„ 

A(t) = 5 i 2 i 

<sin^9>^ - <sia20> ^ 

<S(t,X)sia2e>^ - <S(t,X)>j^<sin2e>2j 

B(t) = —7 (10b) 

<sia^e>^ - <sia2e> ^ 

The angle. 5, in equations 10 is ideally the incident reflection angle at each interface. To convert the offset 
and travel time of the field data to an actual incident reflection angle requires information about interval 
velocities, which may not always be known. However, various approximate solutions have proven to be 
satisfactory. For example, the same functional form, A + Bsin^e holds to order sin^e. whether the angle 
used is the incident reflection angle or some other angle, 0 . as long as the two angles are related by an 
expression of the form 
sm$ = Esina + Fsin^e (n) 

Thus. If an injection angle is employed instead of an incident reflection angle, this expression is still exact in 
the lowest order of d because of Snell's law. Even an average, or effective reflection angle is related to the 
actual reflection angle by an expression of the form of equation 11. 
In practice an effective angle 6 as defined by 



may be used for measurement of the offset dependence of reflection amplitude. In this case Vg (t) Is a 
reference effective velocity* function which may be refined in an iterative manner during the process of 
velocity analysis, as described more fully below. 

If a more detailed velocity function is available it may be used to relate offset and travel time to 
reflection angle for use in equations 10. In this case ray tracing through a velocity model may be 
performed. The methods of the present Invention may employ any type of ray-tracing technique. Further, 
although much of the detailed discussion given herein is couched in terms of velocity models which vary 
only with depth, the inclusion of lateral velocity variations and dipping reflections in the model through 
which the ray-tracing occurs is within the scope of the present Invention. 

From the ray tracing, various values of and sin (one for each two-way travel time and each offset) may 
be calculated and stored as a function of the source-receiver offset distance X. and the two-way travel time 
to the reflector. For purposes of computational economy these values may be calculated for a limited or 
preselected number of offset and time tie points. The values at these tie points may be interpolated to 
determine sine for other offsets and times. An amplitude correction term may also be calculated and stored 
with the foregoing stored data; this amplitude term represents the calculated change in amplitude due to 
geometrical spreading of a wave-front for the various raypaths modeled using ray tracing. Injection and 
emergence angles may also be calculated during the ray tracing and stored for use in additional amplitude 
corrections, due for instance to array effects and absorption. 

The quality of fit coefficient, E. and the normal incidence reflection amplitude estimate. A(t), are both 
insensitive to the manner in which sin^e is estimated. 

The measured offset dependence, B(t), is much more sensitive to this choice but may be corrected at 
the time of interpretation by modelling with a more accurate interval velocity function measured with well 
logs or estimated from an appropriate inversion scheme. 

An advantageous method of the present invention for performing improved moveout velocity analysis 
Initially selects an earth velocity model representative of the area of interest and then employs ray tracing 
through this model to determine sine for a plurality of values of offset and travel time. The data may then be 
con-ected for any significant change in amplitude or in shape due to field instrumentation and travel through 
the model. Examples of such phenomena, are seismometer array effects and Inelastic absorption in the 
earth. The data are also corrected for wavefront spreading amplitude changes as determined during the ray 
tracing, or by an appropriate, separate ray tracing step. 



10 



EP 0 335 450 A2 



Using equation 4, the arrival times of tlie data are then dynamically corrected for each of several trial 
effective velocity functions and the terms of equations 10 calculated. A(t), B(t), and E(t) are computed using 
equations 10 and 7, and tables and/or displays of the quality of fit coefficients. E, for each trial effective 
velocity are made. The effective velocities at the local maxima of the quality of fit coefficients (see for 

5 example Fig. 5A) are called "selected velocities" and are chosen as the best or optimum estimate of 
effective velocities for the times of the maxima, where the estimation process includes changes of signal 
amplitude and shape with offset As with conventional velocity analysis, only one maximum is chosen at any 
value of t and for some values of t no suitable maximum value may be found, which may require 
interpolation between selected maxima as described later herein. 

70 Additional outputs of this method are the least-squared-error fitting coefficients A(t) and B(t) for such 
selected velocities. In the sense that a pair of coefficients A(t) and B(t) are calculated for each trial effective 
velocity function and are determined by the velocity function, the selection of an "optimal" velocity function 
determines the "optimal" set of coefficients A{t) and B(t): this determination is defined herein as the 
"simultaneous" determination of the optimal velocity function and the corresponding optimal Aft) and B(t) 



Although the methods of the present invention Involve calculations in addition to those used in 
conventional velocity analysis, these methods are not dramatically more costly than conventional velocity 
analysis and yield A(t) and B(t) traces in addition to improved velocities. However, for determination of the 
amplitude variation with offset to yield meaningful elastic properties and densities, care must be exercised 
20 during trace amplitude adjustment of the data, e.g., avoiding arbitrary gain correction such as automatic 
gain control (AGC) and making other amplitude corrections as described above. 

Optionally, the methods of the present Invention may be iterated by using the optimum selected 
velocities from one iteration as the starting velocity model in a subsequent iteration of the processing steps 
which will then yield refinements to the earth velocity model and improved estimates of amplitude 
25 adjustments, angles (fl), and optimum Ve (t), A(t), and B(t) values. 

The steps for implementation of an advantageous method of tine present invention may be summarized 



(1) A reference effective velocity function is used to calculate the relation between sm(e) and X and t 
either In an approximate form using equation 12 or by ray tracing based on the corresponding interval 
30 velocity function. The methods for determining average interval velocities from an effective velocity function 
are well known. For computational economy these calculations may be performed for a limited number of 
offset and time tie-points. Later, as needed, these tie-point values may be interpolated to determine sin(fi) 
for other offsets and times. 



(2) Seismic data for a CMP gather Is preprocessed, If required (frequency filtering, static corrections, 



35 amplitude conrections, common offset migration, etc.). Note tiiat although the amplitude correction which is 
commonly applied in conventional moveout velocity analysis is a short gate automatic gain control (AGC) 
designed to make the data fit as well as possible the conventional model assumptions (constant wavelet 
amplitude for any offset), the amplitude adjustment applied for tiie methods of the present invention should 
not involve AGC but includes a wavefront spreading correction as determined by ray tracing using the same 

40 velocity function as described in step (1), as well as corrections for known or estimated channel 
sensitivities, source strengths and receiver coupling variations. Corrections for array effects and absorption 
may be determined using tiie ray tracing results and should also be included in the amplitude adjustment If 
they are significant. 



(3) For each trial effective velocity function tiie data is dynamically corrected in the same manner as 



45 in conventional moveout analysis. For a conventional moveout analysis two sums are computed for each 
time sample: 



75 traces. 



as follows: 



50 



— Z S(t,X) , 



and 



(14a) 



1 



55 



— I S(t,X)2. 



(14b) 
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For the method of the present invention three additional sums are computed: 



1 



(14C) 



I 



and 



X 



1 



I 



sia*(8(X)) , 



which depend upon the offeets (i.e. acquisition geometry) used, but not the data, and 



1 



(146) 



S(t,X) sin' 




which does depend upon the data. 

(4) Finally the least-squared-enror fitting coefficients A(t). B(t) and the quality of fit coefficient E(t) are 
computed using the sums (14a-14e) in equations 7, 9. and 10. Tables and/or displays of the quality of fit 
coefficients are made for each velocity function and time sample and then local maxima of these fit 
coefficients may be picked as in conventional velocity analysis, A velocity function may then be selected to 
pass through these picked maxima, and is referred to as the "optimal" velocity function, 

(5) Additional results of this method are the least-squared-error fitting coefficient traces A(t) and B(t). 
For each trial velocity function one A(t) and one B(t) trace are determined. These traces may be used to 
interpolate optimal A(t) and B(t) traces after the optimal velocity function has been determined, or optimal A- 
(t) and B(t) traces may be recalculated from the optimal velocity function. 

(6) The optimal A(t) and B(t) traces may then be considered, either individually or in combination, as 
characteristic attributes of the earth and these attributes may be displayed and interpreted by conventional 
or other methods. For example, the traces provide information about elastic properties and densities which 
may be interpreted in terms of Ilthology. fluid content, or other geologically significant features. 

The basic steps of the method for processing gathered data to estimate a moveout velocity function, a 
norma! incidence trace, and an offset dependence trace may be briefly summarized by referring to Fig. 3. 
As noted in Fig. 3. there may be previous processing steps before initiation of the method of the present 
invention (block 301). preferably including a gather based on common source-receiver midpoint or common 
reflection point, and appropriate amplitude corrections for known effects. The method initially selects an 
earth model that includes a velocity model based on known or estimated parameters (block 302). This initial 
velocity model is referred to as a "reference" velocity model, Raytracing Is performed through the reference 
model to determine e and sin^ e for various offsets and travel times (block 303). The seismic reflection data 
from the area of interest Is dynamically corrected for each of a plurality of trial moveout velocity functions 
(block 304). Fitting coefficient traces A(t) and B(t) are determined for each zero offset travel time by least- 
squared-error fitting the seismic amplitudes in the data at that time to a form, which is preferably A + Bsin^e 
(block 304). A quality of fit coefficient is calculated (block 304), An "optimal" moveout velocity function is 
selected that interpolates between the velocities corresponding to maxima of the quality of fit coefficient 
(block 305). The optimal normal incidence trace A(t) and offset dependence trace B(t) may then be 
determined by interpolation or by repeating the dynamic correction process with the optimal moveout 
velocity function. 

This process may be repeated for each CMP location along a seismic line to provide a display of A(t) 
and B(t) traces, or combinations thereof, at a plurality of CMP locations, as shown in Rgures 9, and 
discussed later herein. Frequently, for computation economy, the optimal moveout velocity function may not 
be determined at every CMP point and for such a case dynamic correction and A(t) and B(t) trace 
estimation at intermediate CMP locations is accomplished with a moveout velocity function interpolated 
between the optimal functions determined at adjacent CMP locations. 

Optionally, subsequent processing steps may or may not be performed (block 307); these steps may 
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include several iterations of the foregoing process, using the optinnal moveout velocity selected in the latest 
iteration to determine an improved reference model for the next iteration. Thus, it may be seen that the 
present invention provides a method for processing seismic data that simultaneously determines by an 
optimization process the moveout velocity, nomnal incidence reflection amplitude, and angular dependence 

5 of reflection amplitude from the seismic data. 

The optimal A(t) and B(t) traces may then be considered, either individually or in combination; as 
characteristic attributes of the earth and these attributes may be displayed and interpreted by conventional 
or other methods (block 308). For example, the traces provide information about elastic properties and 
densities which may be interpreted in terms of lithology, fluid content, or other geologically significant 

10 features. 

To illustrate the effectiveness of the methods of the present Invention a synthetic data example will be 
considered. In the example, synthetic seismogram traces for various values of X were computed for an 
earth model in which velocity and density vary only with Z. The computation used realistic estimates of 
vertical variation In Vp , Vs , and density taken from actual well logs and a wavelet estimated from 

15 conventional seismic data. Angles of incidence for each time and X were computed by ray tracing and the 
reflection coefficient for each time and X was computed from the Zoeppritz equation for the appropriate 
angle of incidence. The wavelet and the reflection coefficients were then convolved to determine the 
synthetic trace at each time and for each X. 

Fig. 4 depicts this synthetic data dynamically corrected with the true effective velocity computed from 

20 the model. The vertical axis represents time (in seconds), whereas the horizontal axis represents offset (in 
feet). It can be seen, for example, between two-way travel times of 0.9 to 1 .2 seconds in Fig. 4, that the 
conventional approximation that apparent waveform does not vary with offset is not always valid. As was 
indicated earlier In the discussion of Figs. 2A-D. this apparent waveform variation may occur when closely 
spaced reflectors have reflection coefficients with different angular dependences. The effects of the failure 

25 of the assumptions of conventional velocity analysis can be seen in Figs. 5A-B. which depict a comparison 
of conventional velocity analysis and the methods of the present Invention applied to the data of Rg. 4. 

Conventional velocity analysis was performed in the manner described hereinbefore and the results are 
shown in Fig. 5A. Quality-of-fit coefficients were posted at points on a two dimensional grid with trial 
effective velocity and vertical (zero offset) two-way travel time (t) as the two coordinate axes of the grid. For 

30 the displays of these quality of fit coefficients contouring and grey scale shadings were employed. However, 
other displays may be employed. The shading of the grey scale in Figs. 5A-B and 6A-B is representative of 
the magnitude of the fitting coefficients (with a darker shading indicating a higher amplitude than an 
adjacent lighter shading except where the grey scale repeats, I.e.. the light to dark scale repeats to allow 
coverage of a large range of amplitudes of the quality of fit coefficients with a small number of different 

35 shadings). The display techniques used herein are similar to those used in conventional velocity analysis 
but the underlying calculations of the present invention are not the same as those of such conventional 
velocity analysis. 

Referring now to Figs. 5A-B wherein the vertical axes represent time (in seconds) and the horizontal 
axes represent stacking valocity (in ft/sec), there may be seen a comparison of panels of the quality of fit 

40 coefficients from conventional velocity analysis (Rg. 5A) and from the advantageous method of the present 
invention (Fig. SB) for the data of Fig. 4. It may be seen that the maximum quality of fit coefficients for the 
method of the present invention. Fig. 5B, are much more continuous with time than those of the 
conventional velocity analysis, Rg. 5A. Even more importantly, the maxima occur for velocities more closely 
approximating the true effective velocity function determined from the model, which is depicted as curve 

45 "a" on both Rgures. 

Rgs. 6 A and B depict, in a modified display, corresponding portions of the quality of fit coefficient data 
depicted In Rgs. 5A-B, More particularly, Figs. 6 A and B depict panels of the quality of fit coefficients 
recalculated using trial velocities which differ from the true effective velocity by fixed amounts. They are 
plotted against t (vertical axis) and delta Ve, the difference between estimated Ve and true Ve (horizontal 
50 axis). Rgs. 6 A and B show, just below one second, that conventional velocity analysis can result in Ve 
estimate errors as great as t200 ft/sec at times where the method of the present invention gave much more 
accurate values. 

Rgs. 7 A and B depict panels of fitting coefficient traces A(t) and B(t). respectively, for the trial velocity 
functions used to generate Rg. 6B. The vertical and horizontal axes are the same as in Rgs. 6A and B. In 
55 particular, the A(t) and B(t) trace panels. Rgs. 7A and 7B. respectively, consist of A(t) and B(t) traces 
computed with estimates of Ve which differ from the true V^ by the amounts indicated on the delta Ve axis. 
The A(t) and B(t) traces computed with the true Vo are in the center of the panels. The A(t) trace panel. Rg. 
7A. which contains the nonmal Incidence traces, shows that the A traces are relatively Insensitive to velocity 
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errors, that Is, A(t) changes very slowly with delta Ve . From Fig. 7B, it is seen that the B traces are much 
more sensitive to velocity errors and an error of as little as 100 ft/sec may give a poor indication of the true 
offset dependence of the data. Thus, there is a need to select accurate velocities if the offset dependence 
is to be properly characterized. 

5 In Rg. 8 an expanded scale of the B(t) panel, is compared to K. where K is the coefficient of sm^{9) 
computed exactly from the input model and equation 6. (The exact trace Is mari<ed K.) The vertical axis 
represents time (in seconds) whereas the horizontal axis represents delta Vg (feet/sec). There are several 
effects which would cause the measurement of B(t) of the present Invention not to agree precisely with the 
coefficient. K, in equation 6, even when the correct moveout velocity is used, Rrst. the measurement of B(t) 

10 is made with an effective incidence angle derived from a smooth interval velocity function. Second, the 
synthetic seismograms for this data set were generated using Zoeppritz's expression for the offset 
dependence of the reflection coefficients rather than the approximation of equation 6. Third, signal distortion 
due to normal moveout dynamic correction will somewhat limit the accuracy with which offset dependence 
can be measured by this method- (Note that all three of these effects may be accounted for at the 

15 interpretation stage.) In spite of these effects the agreement between the measured B(t) trace with small Ve 
errors and the theoretical K calculated for this earth model is quite good and much better on average than 
the results obtained using Vg from conventional velocity estimation procedures. 

The Interpretation of the normal incidence A(t) trace and offset dependence trace B(t) may be 
performed in many ways. For example, lateral variations in the A(t) and B(t) traces themselves may be 

20 empirically associated with local variation of lithology or fluid content. Other interpretations may best be 
made from various A(t) and B(t) trace combinations and ratios which tend to enhance or discriminate in 
favor of a particular physical property. For instance, in sedimentary rocks at depths of typical exploration 
interest, the Vp/Vs ratio is often near 2. If this value is incorporated into the approximation given by 
equations 6 and 2, one obtains, 

25 Rp(9) = (i=la + l^p) + (Ra-Rp-2R;g)sin2(e) (15a) 

If the small contrast expressions of equation 2f and equation 2g are used then equation 15a further 
simplifies to. Rp(0) = Rp(0) + [Rp(O)-2Rs(O)]sin20. {15b) 

For the approximation of equation 15b, (A-B)/2 should be a measure of the shear-wave reflectivity function 
Rs (O). Lateral changes in compressional-wave reflectivity. Rp , as measured by A, which are not 

30 accompanied by a change in shear wave reflectivity. Rs.. as measured by (A-B)/2, may be a good 
indication of a change in pore fluid. An example of this type of display is shown in Rgs, 9A and 9B. 

Fig. 9A depicts a panel of the A traces for field seismic data (rather than synthetic) at a series of 
adjacent common midpoints and Rg. 9B depicts a panel of corresponding A-B traces for these same 
locations. Both Rgs. 9A and 9B wherein the vertical axes represent time (in seconds) have two seismic 

36 events marked "b" and "h" and common midpoint locations marked Xi, Xa, X3, and X4. Event "b" is a 
consistent right loop on both Figs. 9A and SB, Event "h" has a consistent strong left loop on Fig. 9A 
between X3 and X4 but changes to a weak and variable loop between Xi and Xa. This systematic difference 
(i.e. increase in amplitude) between X3 and X^ Is not visible on Rg. 9B. Thus, it may be suspected that 
event "b" is a reflection from a layer which does not have a difference in fluid content along the line of the 

40 section shown, e.g., it may be consistently brine-filled, while event "h" may be a reflection from a layer in 
which the fluid content differs along the line, in this case probably indicating gas between X3 and X4. 

An indicator of small differences in the Vp A/s ratio on two sides of an interface may be determined by 
making the same approximations as those of equation 15b, along with the assumption that Rp is small. 
Taking the derivative of Vp A/s yields 

45 

Thus, the sum of A + B may be a good indicator of Vp /Vs ratio and lateral changes in the sum may indicate 
changes in rock type or fluid content. In typical cases, the change from brine-filled to gas-filled rocks 
55 produces a significant change in Vp Ns . Further, if the compressional velocities are nearly constant, one 
would expect a change in Vp A/s to be an indication of a change in lithology. 

In both of these cases simplifying approximations have been made which allow one to qualitatively 
predict the cause of observed changes in various displays. In most cases it is desirable to confirm these 
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predictions using more accurate methods. One such method uses standard inverse modelling procedures 
similar to those suggested by Backus and Gilbert. See Backus, G.E. and Gilbert. F.J.. "Numerical 
applications of a formalism for geophysical Inverse Problems". Geophys. J.R. Astr. Soc. 13. pp. 247-276, 
(1967); Backus, G. and Gilbert, F. "The resolving power of gross earth data". Geophys. J.R. Astr. Soc. 16, 
pp. 169-205. (1968): and Backus. G. and Gilbert. P.. "Uniqueness in the inversion of Inaccurate gross earth 
data", Phil, Trans. Roy. Soc. London 266. pp. 123-192, (1970). 

Hence from the above, it is clear that the present invention provides methods for measuring effective 
velocity while appropriately accounting for the varying offset dependence of subsurface reflectors. The 
present invention determines traces A(t) (the "normal incidence trace") and B(t) (the offset dependence of 
the seismic data) which are desirable for useful lithologic and/or fluid content interpretation of seismic data. 

The specific sequence of steps described hereinbefore for the methods of the present Invention may be 
changed and stlH achieve the same result without departing from the scope of the present invention. As 
noted hereinbefore, the methods of the present invention are applicable to both compressional wave and/or 
shear wave seismic exploration techniques. 

Many other variations and modifications may be made in the techniques hereinbefore described, by 
those having experience in this technology, without departing from the concepts of the present invention. 

Accordingly, it should be clearly understood that the methods depicted in the accompanying drawings 
and referred to in the foregoing description are illustrative only and are not intended as limitations on the 
scope of the invention. 

Claims 

1. A method for processing seismic data to enhance information In seismic signals reflected from 
contrasts or differences In elastic constants or densities in the subsurface of the earth, comprising the steps 
of: 

- generating seismic energy on or near the surface of the earth; 

- measuring on or near said surface seismic signals resulting from said seismic energy reflected from 
interfaces, and recording said signals to obtain seismic data characterized by the step of specifically 
including the effects of amplitude variations with offset, 

2. The method as claimed in claim 1. characterized by simultaneously determining moveout velocity, 
normal incidence reflection amplitude, and offset dependence of reflection amplitude from said data. 

3. The method as claimed in claim 2. characterized in tiiat said simultaneously determining from said 
data moveout velocity, normal incidence reflection amplitude, and offset dependence of reflection amplitude 
takes place by an optimization process which recognizes the interdependence of said moveout velocity and 
offset dependence of reflection amplitude. 

4. The method as claimed in claim 1 . characterized by: 

- calculating from a known reference velocity function effective reflection angles for various two-way travel 
times and offsets of said data, 

- determining from said data and effective reflection angles, for a trial moveout velocity function, fitting 
coefficients for a preselected offset dependence function which expressly accounts for offset dependence of 
reflection amplitudes, 

- repeating said determining of said fitting coefficients for a plurality of said trial moveout velocity functions, 
and 

- selecting for various normal incidence two-way travel times appropriate moveout velocities and fitting 
coefficients. 

5. The method as described in claim 4. characterized by determining seismic attributes from said 
selected fitting coefficients. 

6. The method as described in claim 4, characterized by determining an improved reference velocity 
function from said selected appropriate moveout velocities, and redetermining said fitting coefficients and 
moveout velocities employing said improved reference velocity function. 

7. The method as described in claim 6. characterized by displaying for said selected moveout velocities 
said fitting coefficients that are representative of measured normal incidence and measured offset depen- 
dence for said data. 
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8. The method as claimed in claim 1 , characterized by determining from said measured seismic signals 
in a manner which includes amplitude variation with offset optimized moveout velocities, normal incidence 
reflection amplitudes, and offset dependence of reflection amplitudes, determining elastic properties and 
densities of said subsurface materials from said moveout velocities, normal Incidence reflection amplitudes. 

5 and offset dependencies of reflection amplitudes. 

9. The method as claimed in claim 1 , characterized by the steps of: 

- selecting a reference velocity function for a region of interest, 

- calculating reflection angles for various acquisition offsets and two-way travel times used to acquire said 
seismic data. 

70 - determining moveout velocities by optimizing a function expressly including offset dependence of 
reflection amplitudes, 

- for the determined moveout velocities simultaneously determining coefficients that fit the offset depen- 
dence of said data, and 

- displaying for said selected moveout velocities attributes determined from said coefficients that are 
75 representative of measured zero offset and measured offset dependence for the seismic data. 

10. The method as claimed in claim 1. characterized by the step of determining moveout velocities in a 
manner that includes amplitude variation with offset. 

11. The method as claimed In claim 1, characterized by the step of determining normal Incidence 
reflection amplitude in a manner that includes amplitude variation with offset. 

20 12. The method as claimed in claim 1, characterized by the step of determining amplitude variation with 
offset in a manner that includes simultaneous moveout velocity analysis. 
13. The method as claimed in claim 1, characterized by the steps of: 

- developing a velocity model having known or estimated parameters, 

- performing raytracing through said model to detenmlne Q and sin^e for a plurality of offsets and travel 
25 times. 

- dynamically correcting said data with a plurality of selected moveout velocity functions, 

- determining fitting coefficient traces corresponding to normal Incidence reflection amplitude and amplitude 
variation with offset for each zero offset travel time of said data by fitting the amplitudes of said seismic 
data to a form A + Bsin^e using a least-squared-error fitting, 

30 - calculating a quality of fit coefficient for said least-squared-error fitting, 

- selecting a moveout velocity function that maximizes the quality of fit coefficient for a plurality of travel 
times. 

- interpolating a continuous moveout velocity function between said plurality of travel times, and then 

- determining fitting coefficient traces corresponding to said interpolated moveout velocity function. 
35 14. The method as claimed in claim 13. characterized by the steps of: 

- modifying said velocity model based upon said selected moveout velocity function, and then 

- reperforming the foregoing steps to select an Improved moveout velocity function, and then 

- redetermining fitting coefficient traces corresponding to said improved moveout velocity function. 

15. The method as claimed In claim 14. characterized by the step of displaying said fitting coefficient 
40 traces. 

16. The method as claimed in claim 14, characterized by the step of determining seismic attributes 
from said fitting coefficient traces. 

17. The method as claimed in claim 16, characterized by the step of displaying said seismic attributes. 

18. The method as claimed in claim 16, characterized by the step of interpreting s^d seismic attributes 
45 to detenmlne elastic properties and densities. 
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